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^ ■ Abstract: Computational modehng and simulation are presented on the motion of red blood 

cells behind a moving interface in a capillary. The methodology is based on an immersed 
boundary method and the skeleton structure of the red blood cell (RBC) membrane is 
i__il modeled as a spring network. The computational domain is moving with either a designated 

Ch ' RBC or an interface in an infinitely long two-dimensional channel with an undisturbed flow 

i-rt*' field in front of the domain. The tanking-treading and the inclination angle of a cell in a 

i . simple shear flow are briefly discussed for the validation purpose. We then present the results 

r— I ! of the motion of red blood cells behind a moving interface in a capillary, which show that 

j/5 ! the RBCs with higher velocity than the interface speed form a concentrated slug behind the 

y ' interface. 

C/3 . 
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> 1 Introduction 

00 ■ 

(N; 

?^ ■ The rheological property of the red blood cells (RBCs) is a key factor of the blood flow char- 

• ■ acteristics at the microchannel level, especially the particulate nature of the blood becomes 

Q ! significant when studying blood drop through a glass capillary within miniature blood diag- 

^^ I nostic kit. The penetration of the blood suspension in a perfectly wettable capillary has been 

• • I analyzed in [H [2] . The failure of such penetration is attributed to three RBCs segregation 

.^ ■ mechanisms: (i) corner deflection at the entrance, (ii) the intermediate deformation- induced 

r> ■ radial migration and (iii) shear-induced diffusion within a packed slug at the meniscus. The 

c^ , key mechanism responsible for penetration failure is the deformation-induced radial migra- 

tion, which endows the blood cells with a higher velocity than the meniscus to form the 
concentrated slug behind the meniscus (see Figured]). The results in [1] |2] shed light on 
making the smallest microfluidic kit and loading microneedle that require the least amount 
of blood sample. 

Nowadays in silico mathematical modeling and numerical study of RBC rheology have 
attracted growing interest (see, e.g., [21 H]). The immersed boundary method developed by 
Peskin, e.g, [5l[71l6], has been one of the popular methodologies for numerically studying the 
RBC rheology due its distinguish features in dealing with the problem of fluid flow inter- 
acting with a flexible fluid/structure interface. For example, in [S]-[T7], immersed boundary 
methods have been combined with different RBC membrane models to simulate the motion 
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Figure 1: Schematics of the BRCs moving behind a meniscus. 



of RBCs and vesicles in fluid flow. In [151 [THl EI], we have successfully combined an im- 
mersed boundary method with a spring model developed in [18] to simulate the motion of 
RBCs in shear flows and Poiseuille flows. In this paper we have generalized the aforemen- 
tioned methodology to simulate the RBCs aggregation behind a moving interface considered 
in [H [2] by having the computational domain moving with an interface in an infinitely long 
two-dimensional channel with an undisturbed flow field in front of the domain since the 
typical periodic boundary condition in the direction of the channel wall is not well suited 
anymore. To mimic the motion of the RBCs behind a meniscus in a capillary, we have 
considered a flat interface moving with a given constant speed in this paper. The simulating 
results of the motion of red blood cells behind a moving interface show that the RBCs with 
higher velocity than the interface speed form the concentrated slug behind the interface, 
which resembles the motion of the RBCs observed in [H [2|. The structure of this paper is as 
follows: We discuss the elastic spring model and numerical methods in Section 2. In Section 
3, the tanking-treading and the inclination angle of a cell in a simple shear flow are briefly 
discussed for the validation purpose. We then present the results of the motion of red blood 
cells behind a moving interface in a capillary. The conclusions are summarized in Section 4. 



2 Models and methods 

Let f2 be a bounded rectangular domain filled with blood plasma which is incompressible, 
Newtonian, and contains RBCs with the viscosity of the cytoplasm same as that of the blood 
plasma (see Figure [2]). For some T > 0, the governing equations for the fluid-cell system are 
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V-u = in f], t G (0,T) 



(2.2) 



where u and p are the fluid velocity and pressure, respectively, p is the fluid density, and /i 
is the fluid viscosity, which is assumed to be constant for the entire computational domain. 




Figure 2: An example of the computational domain with a cell. 



In ( 12. ip . f is a body force which accounts for the force acting on the fluid/cell interface. 
Equations (12. ip and (12. 2p are completed by the following boundary and initial conditions: 



u = go on Yd, 

du 

11- np = on r„, 

an 

u(0) = uo 



(2.3) 
(2.4) 
(2.5) 



where the domain Vt is taken from an infinitely long channel with its boundary denoted by 



r = utiF, 

r„ = and Td = T, (ii) r„ = T^ and Td -- 
Poiseuille flow or simple shear flow on Td 



In the simulations, we have considered two types of boundary conditions: (i) 

Fi U r2 U F3 with go having the profile of either 



2.1 Elastic spring model for the RBC membrane 



A two-dimensional elastic spring model used in [18] is considered in this paper to describe the 
deformable behavior of the RBCs. Based on this model, the RBC membrane can be viewed 
as membrane particles connecting with the neighboring membrane particles by springs, as 
shown in Figure [31 Elastic energy stores in the spring due to the change of the length / of 
the spring with respected to its reference length Iq and the change in angle 6 between two 
neighboring springs. The total elastic energy of the RBC membrane, E = Ei + E^, is the sum 
of the total elastic energy for stretch/compression and the total energy for bending which, 
in particular, are 

Ei = ^y,{^-^r (2.6) 



and 



i=l 



N 



k 



E, = ^J2tan\e,/2) 



(2.7) 



In equations (12. 6 p and (12. 7p . A^ is the total number of the spring elements, and ki and kb 
are spring constants for changes in length and bending angle, respectively. 




Figure 3: The elastic spring model of the RBC membrane 

Remark 2.1. In the process of creating the initial shape of RBCs described in [18], the RBC 
is assumed to be a circle of radius i?o = 2.8 fim initially. The circle is discretized into N = 7Q 
membrane particles so that 76 springs are formed by connecting the neighboring particles. 
The shape change is stimulated by reducing the total area of the circle through a penalty 
function 

r. = ^C-^r (2.8) 

where s and s^ are the time dependent area of the RBC and the equilibrium area of the 
RBC, respectively, and the total energy is modified as E + Tg. Based on the principle of 
virtual work the force acting on the ith membrane particle now is 

F. ^ -?^±^ (2.9) 

where r^ is the position of the ith membrane particle. When the area is reduced, each RBC 
membrane particle moves on the basis of the following equation of motion: 

mf, + 7r, = Fi (2.10) 

Here, () denotes the time derivative; m and 7 represent the membrane particle mass and the 
membrane viscosity of the RBC. The position Fj of the ith membrane particle is solved by 
discretizing f l2.10p via a second order finite difference method. The total energy stored in 
the membrane decreases as the time elapses. The final shape of the RBC is obtained as the 
total elastic energy is minimized (please see [I9]). The area of the final shape has less than 
0.001% difference from the given equilibrium area Se and the length of the perimeter of the 
final shape has less than 0.005% difference from the circumference of the initial circle. The 
reduced area of a RBC in this paper is defined by s* = Se/TiRl- 

Remark 2.2. When simulating the case involving a moving interface, we have applied a 
repulsive force to prevent the overlapping between cell and wall. The repulsive force is 
obtained from the following Morse potential (e.g., see [20] ) 



4>{d) = kr{l - e-('^-°'«))2 



where the parameter d is the shortest distance between the membrane particle and the wall 
and do is the range of the repulsive force (when the distance d is greater than do, there is no 
repulsive force). The parameter k,. is a constant for the strength of the potential. 

2.2 Immersed boundary method 

The immersed boundary method developed by Peskin, e.g, [51 dim, is employed in this study 
because of its distinguish features in dealing with the problem of fluid flow interacting with a 
flexible fluid/structure interface. Over the years, it has demonstrated its capability in study 
of computational fluid dynamics including blood flow. Based on the method, the boundary 
of the deformable structure is discretized spatially into a set of boundary nodes. The force 
located at the immersed boundary node Tj = (rji,rj2) affects the nearby fluid mesh nodes 
X = (xi,X2) through a 2D discrete Wunction /^^(x — Fj): 

f(x) = ^F,Z}/,(x-r,) /or |x-r,|<2/i, (2.11) 

where h is the uniform finite element mesh size and 

^/.(x - Ti) = 6h{xi - ri^i)6h{x2 - ri,2) (2.12) 

with the ID discrete 5-functions being 

■^ ^3 - 2\z\/h + ^l + A\z\/h-A{\z\/hy^ , 1^1 < h, 

Sh{z) = <; ^ (5 - 2\z\/h - ^-7 + 12\z\/h-4{\z\/hy) , h < \z\ < 2h, (2-13) 
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0, otherwise. 

The velocity of the immersed boundary node Fj is also affected by the surrounding 
fluid and therefore is enforced by summing the velocities at the nearby fluid mesh nodes x 
weighted by the same discrete (5-function: 

U(r,) = ^/iVx)D,,(x-r,) /or |x-x,| <2/i. (2.14) 

After each time step, the position of the immersed boundary node is updated by 

rri = rr + AtU(rr). (2.15) 

2.3 Space approximation and time discretization 

Concerning the finite element based space approximation of {u,p} in problem (I2.ip - (l2.5p . 
we use the Pi-iso-P2 and Pi finite element approximation (e.g., see [22] (Chapter 5)). For a 
rectangular computational domain Q G R"^ , let Th he a finite element triangulation of Q for 
velocity and 72h a twice coarser triangulation for pressure where /i is a space discretization 
step. We introduce the finite dimensional spaces: 



Woh = {v/i|v/i G Wh, Vft, = on r^}, 

Ll = {qhlQh e C\U),qh\T G Pi,VT G Ts, }, 

{g/ilgh G L^, / qhdyi = 0} 



7^2 



where Pi is the space of polynomials in two variables of degree < 1. We apply the Lie's 
scheme [211 [22] with the above finite elements to equations fl2.ip - fl2.5p with the backward 
Euler method in time for some subproblems and obtain the following sequence of fractional 
step subproblems (some of the subscripts h have been dropped): 

u° = Uq is given; for n > 0, u" being known, we compute the approximate solution via 
the following fractional steps: 



1. Update the position of the membrane by fl2.14p and fl2.15p and then compute the force 
P based on the fluid/cell interface by (ESI) and (l2lll) . 



2. Solve 



/ ^^ ■ vdx + / (u" ■ V)u(t) ■ vdx = 0, on (t", t 

u(r) = u", 

u(t) G Wh, u(t) = go,,, on r- X (r,r+i), 



n j.n+1^ 



VvG w; 



and set u"+2/3 = u(t"+^). 
3. Finally solve 

U"+1 _ ^"+2/3 
At 



vrfx + /i / Vu''+^ : Vvdx 



Jn Jn 



gV-u"+Mx = 0, VgGL 



h' 



Oft,' 



(2.16) 



(2.17) 



I u"+i G W^,, u(t) = go,,, on r,; p"+i G Ll (p"+i G L^ ^ if r„ = 0). 



In eq. 02.161). we have T^ = {x|x G T, go,h(x) ■ n(x) < } and H^gh = {^hWh e lVh,v,j = 
on r^}. The quasi- Stokes problem (12.171) is solved by a preconditioned conjugate gradient 
method (see, e.g., [22] )• The subproblem (I2.16p is an advection type subproblem. It is solved 
by a wave-like equation method, which is described in detail in [23] and ^^. 

Remark 2.3. In simulations, the computational domain Q moves to the right with either 
the mass center of a RBC or the interface (see, e.g., [2S1 [2S] and references therein for 
adjusting the computational domain according to the position of the particle). Due to the 



use of structured and uniform mesh in our simulations, it is relatively easy to have the 
computational domain moving with a designated cell. Generally when the mass center of a 
RBC moves to the right in an infinitely long channel, we add one vertical grid line to the 
right end of the computational domain if the cell mass center crosses one vertical grid line 
after we predict its new position and at the same time we drop one vertical grid line at the 
left end of the computational domain. In the mean time at these new grid points added 
at the right end, we assign the values of velocity field according to either Poiseuille flow or 
simple shear flow depending on the test cases. When following an interface moving to the 
right with a constant speed, we have applied the same strategy. 



3 Numerical results and discussion 



3.1 Tank-treading of a single cell in shear flow 
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Figure 4: (Color online). Steady inclination angle versus the cell swelling ratio (left) and 
membrane tank-treading velocity (scaled by 7-R0/2) versus the cell swelling ratio (right) in 
comparison to Shi. et al. [15] and Kaoui. et al. [I3] in different cases. Case I: 112/im x 7/im 
domain with Dirichlet boundary conditions. Case II: 80yum x 7/im domain with Dirichlet 
boundary conditions. Case III: 80/im x 7/im domain with Neumann inflow condition and 
Dirichlet outflow condition. 



We have first validated the computational methodology with two types of boundary 
conditions discussed in Section [2] by comparing the inclination angle and the tank-treading 
frequency of a single RBC in shear flow. Here are the parameters used in the simulations: The 
values of parameters for modeling cells are same with [151 [HI [IZ| as follows: The bending 
constant is /cj, = 5 x 10^^*^N ■ m, the spring constant is /;;; = 5 x 10~'^N ■ m, the penalty 
coefficient is k^ = 10~^N ■ m, the repulsive force coefficient is kr = 10~^N ■ m, and the range 
of the repulsive force is do = 2h where h is the mesh size for the flow velocity field. The cells 
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Figure 5: (Color online). The positions of 12 cells in a capillary behind a moving interface 
at t = 0.01, 0.5, 2.5, 4.5 and 6.25 ms and the velocity field with 12 cells at t =6.25 ms (from 
top to bottom). 
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Figure 6: (Color online). The position of 68 cells in a capillary behind a moving interface 
at t = 0.01, 1.8, 5.37 and 10 ms and the velocity field with 68 cells at t = 5.37 and 10 ms 
(from top to bottom). 
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are suspended in blood plasma which has a density p = l.OOg/cm^ and a dynamical viscosity 
/i = 0.012g/(cm ■ s). The viscosity ratio which describes the viscosity contrast of the inner 
and outer fluid of the RBC membrane is fixed at 1.0. The dimensions of the computational 
domain are 112/im x 7yum and 80yum x 7yum. Then the degree of confinement (2Rq/H) is 0.8 
where H is the height of the channel. The grid resolution for the computational domain is 
80 grid points per lO/im. The time step At is 1 x 10~^ms. The initial position of the mass 
center of the cells are (56, 3.5) and (40, 3.5) for the longer domain and the shorter domain, 
respectively. To have a shear flow, a Couette flow driven by two walls at the top and bottom 
which have the same speed U/2 but move in directions opposite to each other is applied to 
the suspension, where the speed U is given hj U = '-f*H with a given shear rate 7. The shear 
rate used in the simulation is 7 = 275/s. The steady inclination angles of the tank-treading 
for four values of s*=0.6, 0.7, 0.8 and 0.9 are presented in Figure IH which show the very 
good agreement with the lattice-Boltzmann simulation results in [13] and those previously 
obtained with periodic boundary conditions in [15]. The membrane tank-treading velocity 
(scaled by ■yRo/2) is also in good agreement with the results in [131 [IS]- The results show 
that there is no significant difference when having the Dirichlet boundary conditions on F 
with the length L = 112 and 80 /im or the conditions (12. 3 p and (12. 4 p on the boundary of the 
shorter domain. 



3.2 Multi-cell aggregation in a capillary behind a moving interface 

For the cases involving a moving interface in a capillary, we have considered the one moving 
to the right with constant speed U to mimic the motion of the RBCs behind a meniscus in 
a capillary. Then the associated boundary condition in (12. 3 p on F^ is go = on Fi U F3 
and go = {U, 0)* on F2 and the boundary condition (12. 4 p is satisfied on F4. We have kept 
all the related parameters the same except the following. We have first considered the case 
of 12 cells of swelling ratio s*=0.481 in a capillary of the height 10/im. The computational 
domain Q is 80/im x 10/im. The interface speed is U = 8/3 cm/s. The cells in the center 
of the channel move faster than those next to the top and bottom walls do due to fact that 
the velocity field behaves like Poiseuille flow as the fluid flow is away from the interface 
and the speed of the interface is slower than the velocity of the fluid flow in the channel 
central region away from the interface (see the velocity field in Figures |5] and [6]) . For the 
cells moving away from the interface, they move back to the central region of the channel 
due to the lateral migration of the cells in a flow field like the Poiseuille flow and then move 
toward the interface. Thus the cells form a slug behind the moving interface and move with 
the interface as in Figure O For the case of 68 cells of swelling ratio s*=0.481 in a capillary 
of the height 20/im, we have considered the computational domain Q = 160/im x 20/im. The 
interface speed is f/ = 8/3 cm/s. These 68 cells behave similarly behind the moving interface 
like the motion of the 12 cells considered in the previous case. But it is much clearly for us to 
see that the cells in the channel central region move faster to the right due to the relatively 
faster flow field. Then the cells are piled up behind the interface and move with the interface 
in Figure O 
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4 Conclusions 

In summary, we have developed computational modeling and methodologies for simulating 
the motion of many RBCs in a capillary behind a moving interface in this paper. The 
methodology is based on an immersed boundary method and the skeleton structure of the 
red blood cell (RBC) membrane is modeled as a spring network. The computational domain 
is moving with either a designated RBC or an interface in an infinitely long two-dimensional 
channel with an undisturbed flow field in front of the domain. The tanking-treading and 
the inclination angle of a cell in a simple shear flow are briefly discussed for the validation 
purpose. The results of the motion of red blood cells behind a moving interface in a capillary 
show that the RBCs with higher velocity than the interface speed form a concentrated slug 
behind the interface, which is consistent with the results in [H |2]. The lateral migration is 
also a key factor for the formation of a slug behind the moving interface. 
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